## Import data from .csv
cont_snow_data <- read.csv("data/snotel_cont.csv")

# SNOTEL Analysis
cont_snow_data_spring <- cont_snow_data %>%
  mutate(date = as.POSIXct(date, format = "%Y-%m-%d"),yday = yday(date), year = year(date)) %>%
  filter(state != "SD", year >= 1979, ifelse(leap_year(year) == FALSE, yday >= 60 & yday <= 243, yday >= 61 & yday <= 244)) %>%
  mutate(newsnow = ifelse(snow_water_equivalent>lag(snow_water_equivalent), snow_water_equivalent-lag(snow_water_equivalent), 0))

# Summarize the new snow fall events by state
snotel_state <- cont_snow_data_spring %>%
  group_by(state, year) %>%
  summarise(newsnow_days = sum(newsnow>0, na.rm = TRUE), n = length(unique(site_id)), newsnow_site = newsnow_days/n)
## `summarise()` has grouped output by 'state'. You can override using the
## `.groups` argument.
# 
# spring_snow_state <- ggplot(snotel_state) +
#   geom_line(aes(x = year, y= newsnow_site, color = state))
# 
# ggplotly(spring_snow_state)

# Summarize the new snow fall events by site
snotel_sites <- cont_snow_data_spring %>%
  group_by(site_id, site_name, year, state) %>%
  summarise(newsnow_days = sum(newsnow>0, na.rm = TRUE))
## `summarise()` has grouped output by 'site_id', 'site_name', 'year'. You can
## override using the `.groups` argument.
# spring_snow_state <- ggplot(snotel_sites) +
#   geom_line(aes(x = year, y= newsnow_days, color = state))
# 
# ggplotly(spring_snow_state)

# Determine average number of days with increased SWE per spring
cont_snotel_site <- cont_snow_data_spring %>%
  group_by(site_id, site_name, state, latitude, longitude, elev) %>%
  summarise(newsnow_days = sum(newsnow>0, na.rm = TRUE), years = length(unique(year)), avg_spring_days = newsnow_days/years)
## `summarise()` has grouped output by 'site_id', 'site_name', 'state',
## 'latitude', 'longitude'. You can override using the `.groups` argument.
mapview(cont_snotel_site, ycol = "latitude", xcol = "longitude", zcol = "avg_spring_days", crs = 4326, grid = FALSE)
cont_snotel_site_yr <- cont_snow_data_spring %>%
  group_by(site_id, site_name, state, year, latitude, longitude, elev) %>%
  summarise(newsnow_days = sum(newsnow>0, na.rm = TRUE))
## `summarise()` has grouped output by 'site_id', 'site_name', 'state', 'year',
## 'latitude', 'longitude'. You can override using the `.groups` argument.
spring_snow_days <- ggplot(cont_snotel_site_yr, aes(x = state, y= newsnow_days, color = state)) +
  geom_boxplot()+
  ggtitle("Days of increased SWE at all sites (1978-2022)")

ggplotly(spring_snow_days)
# Testing things
# test_snow_data <- snotel_test %>%
#   mutate(date = as.POSIXct(date, format = "%Y-%m-%d"),yday = yday(date), year = year(date))
# 
# test_spring <- test_snow_data %>%
#   filter(ifelse(leap_year(year) == FALSE, yday >= 60 & yday <= 243, yday >= 61 & yday <= 244)) %>%
#   mutate(newsnow = ifelse(snow_water_equivalent>lag(snow_water_equivalent), snow_water_equivalent-lag(snow_water_equivalent), 0))
# 
# site_stats <- test_spring %>%
#   group_by(site_id, site_name, state, county, latitude, longitude, year) %>%
#   summarise(day_SN = sum(newsnow>0, na.rm = TRUE))
# 
# ggplot(site_stats) +
#   geom_line(aes(x = year, y= day_SN, color = site_name))